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TECHNICAL MEMORANDUM 


USERS GUIDE FOR SKYLAB DYNAMICS PROGRAM, SKYDYN 

I. INTRODUCTION 


The program SKYDYN was developed to simulate the orbital dynam- 
ics of an uncontrolled asymmetric vehicle subjected to perturbing torques 
due to gravity gradient and aerodynamic forces. The program utilizes 
an oblate rotating Earth model and a variable step size, five-pass Runge- 
Kutta integration scheme. Quaternions are used to represent the attitude 
of the vehicle; thus, there are no restrictions on attitude or small angle 
motion . 

The versatility of the program output allows the user to specify 
the output parameters without reprogramming. Tape output can be used 
for plotting or as input to other programs for open-loop calculations that 
Vtfould otherwise increase the run time for the dynamics program. 

The input data requirements and a sample data input listing are 
given in Section II. Program output specifications and a sample program 
printout are given in Section III. 

Appendix A presents definitions of the program variables. The 
reference frames and corresponding transformation are given in Appendix 
B. The equations of motion used in the simulation are presented in 
Appendix C. A listing of the main routine, SKYDYN, and all subroutines 
used are given in Appendix D. 


II. INPUT DATA REQUIREMENTS 


A . Data Input 

Table 1 defines the variables required as input to the program 
SKYDYN. A change in code number indicates the beginning of a new 
line of input. Each variable specified to be printed and/or saved on tape 
(see Code 2) begins a new line of input with a maximum of 90 lines per 
run. 


All angular data are input in degrees and converted to radians at 
the start of each simulation. A more detailed description of the table 
input for density and aerodynamic coefficients is presented following the 
• list of input data requirements. The transformation matrix from principal 
to body axes (Code 25) is used for open-loop calculations only, and if 
data are unavailable, dummy variables may be used as input. 

All data listed in Table 1 are required for initialization of the pro- 
gram. If it is desired to restart the simulation with initial conditions 
saved from a previous run (IOPTl=l), card input data for Codes 23 and 
24, the initial start time and integration step size are ignored. For 
multiple cases (NCASE > 1, Code 1), data cards Codes 17 through 25 are 
repeated. Card input data are assigned to Unit 5. Tape input (see 
lOPTl, Code 19) is assigned to Unit 14. 

A typical data input listing is given in Section II. D. 


TABLE 1. SKYDYN INPUT 



Program 

Variable 

Units 

Format 

Code 

Symbol 

Definition 

1 

NCASE 

Case number 

Unitless 

15 

(Col 1-5) 


NX 

Number of integration 
variables 

Unitless 

15 

(Col 6-10) 

2 

IVAR 

Location of variable in 
common to be printed and/ 
or saved on tape 

Unitless 

13 

(Col 1-3) 


ISCAL 

Scale factor designation 
ISCAL=0, Scale factor=l.; 
ISCAL=1, Scale factor= 
57.29578; ISCAL^O or H, 
Scale factor supplied by 

Unities s 

11 

(Col 5) 


TABLE 1, (Continued) 


Prog^ram Variable 


Code 

Symbol 

Deflnition 

Units 

Format 


SCAL 

Scale factor 

User 

Specified 

E12.5 
(Col 7-18) 


PNAME 

Name to bo ?*ssigned to 
printed variables 

I’nitless 

A4 

(Col 19-22) 


IPW 

Piint and/or save desig- 
nation: =0, print and 
save; =1, print only; =2 
save only 

Jnitless 

11 

(Col 23) 

3 


Blank card to signify end 
of output variables 



4 

ATNM(I) 

1=1,18 

Atmosphere Title Card 

Unitless 

18A4 

5 

TLAT(I) 

1=1,11 

Latitude table for density 
lookup 

deg 

7F8.3 

6 

TLNG(I) 

1=1,37 

Longitude table for density 
lookup 

deg 

7F8. 3 

7 

FRHO(I,J) 

1=1,11 

J=l,37 

Atmospheric density 

kg/m^ 

6E9.3 

8 

NALP 

Number of total angle-of- 
attack values 

Unitless 

15 

(Col 1-5) 


NPHIA 

Number of aerodynamic 
roll angle values 

Unitless 

15 

(Col 6-10) 

9 

TALP(I) 
1=1, NALP 

Total an gle-of- attack table 
for aerodynamic data 
lookup 

deg 

7F8.3 

10 

TPHIA(I) 
1=1, NPHIA 

Aerodynamic roll angle 
table for aerodynamic 
data lookup 

deg 

7F8.3 

11 

FCA(I,J) 
1=1, NPHIA 
J=1,NALP 

Axial force coefficient 

Unitless 

7F8.3 

12 

FCN(I,J) 
1=1, NPHIA 
J=1,NALP 

Normal force coefficient 

Unitless 

7F8.3 

13 

FCY(I,J) 

Side force coefficient 

Unitless 

7F8.3 


I=1,NPHIA 

J=1,NALP 


3 




TABLE 1. (Continued) 



Code 

Program 

Symbol 

Varioblc 

Definition 

Units 

Format 

14 

FCM(I,J) 

IbI.NPHIA 

J=1,NALP 

Pitching moment coefficient 

Unitless 

7F8.3 

IS 

FCEN(I,J) 

I:=1,NPHIA 

J-1,NALP 

Yawing moment coefficient 

Unitless 

7F8.3 

16 

FCL(I,J) 

I=1,NPHIA 

J=1,NALP 

Rolling moment coefficient 

Unitless 

7F8.3 

17 

CASE(I) 

1=1,18 

Case title card 

Unitless 

18A4 

18 

PROPT 

Multiplier for print inter- 
val (Print interval-PROPT 
♦DTP) 

Unitless 

FlO.l 

19 

DT 

Initial integration time 
step 

sec 

F10.4 


DTP 

Output frequency for 
save tape 

sec 

F10.4 


DTSAM 

Specified time to save 
variables for restart 

sec 

F10.4 


TRUN 

Total run time 

sec 

F10.4 


TIME 

Initial start time 

sec 

F10.4 


lOPTl 

Initialization option 
=0, use initial conditions 
from data pack 
=1, read in initial 

conditions from tape 

Unitless 

15 

20 

WT 

Vehicle weight 

lb 

F10.4 


XCG 

Vehicle eg in x-direction 

ft 

F10.4 


YCG 

Vehicle eg in y-direction 

ft 

F10.4 


ZCG 

Vehicle eg in z-direction 

ft 

F10.4 


XMRP 

Aerodynamic moment refer- 
ence point in x-direction 

ft 

F10.4 


YMRP 

Aerodynamic moment refer- 
ence point in Y-direction 

ft 

F10.4 


ZMRP 

Aerodynamic moment refer- 
ence point in z-direction 

ft 

F10.4 


4 . 


TABLE 1. (Continued) 


Code 

Program 

Symbol 

Variable 

Definition 

Units 

Format 

21 

DREF 

Aerodynamic reference 
diameter 

ft 

F10.4 


SREF 

Aerodynamic reference 
area 

ft^ F10.4 

22 

IXYZ(1,1) 

Moment of inertia about 
x~axis 

slugs-ft^ FlO.l 


IXYZ(2,2) 

Moment of inertia about 
y-axis 

slugs-ft^ FlO.l 


IXYZ(3,3) 

Moment of inertia about 
z-axis 

slugs-ft^ FlO.l 


IXYZ(1,2) 

xy product of inertia 

slugs-ft^ FlO.l 


IXYZ(1,3) 

xz product of inertia 

slugs-ft” FlO.l 


IXYZ(2,3) 

yz product of inertia 

slugs-ft^ FlO.l 

23 

PSI 

Geocentric latitude (posi- 
tive north and negative 
south of equator) 

deg 

F10.6 


LAMDE 

Earth fixed longitude 
(positive east and nega- 
tive west of Greenwich) 

deg 

F10.5 


RMAG 

Radius vector magnitude 

ft 

F10.2 


VIMAG 

Inertial velocity 

ft /sec 

F10.3 


SIGI 

Inertial heading (pos- 
itive clockwise from 
north) 

deg 

F10.6 


GAMI 

Inertial flight path angle 
(positive up from local 
geocentric horizontal) 

deg 

F10.8 

24 

PHIBI 

Initial bank angle 

deg 

F10.5 


ALPHAI 

Initial total angle-of- 
attack 

deg 

F10.5 


PHIAI 

Initial aerodynamic roll 
angle 

deg 

F10.5 


TABLE 1. (Continued) 



Program 

Variable 



Code 

Symbol 

Definition 

Units 

Format 


PQR(I) 

Initial rates about x, y, 

deg/sec 

3F10.5 


1=1 r3 

and z body axes, 
respectively 



25 

ABP(I,J) 

Transformation matrix 

Unitless 

610.5 


J**l,3 from principal axes to 

1=1 , 3 body axes 


26 


Blank card to signify 
termination of input data 


B. Density Table Input 

The density table input is a bivariate function of latitude (tp) and 
longitude (A). The density values were determined for a specified altitude 
using the Jacchia 1970-8 atmosphere model and predicted solar and geom- 
agnetic data provided by Space Sciences Laboratory, Marshall Space Flight 

Center. 

The dependent density values are input in the following manner: 
p(’’Pl> Ap p(\p2> A^) p(ipQ, Xp 

p(^^, A^ 


^37^ ^37^ 


p(i^Q> Ag,^) 


p('l^y> ^37^ 


^37^ 


C 




C. Aerodynamic Coefficients Xable Jfiput 

The aerodynamic coefficients are bivariate functions of aei-odynomic 
roll angle and total angle-of-attack The dimension of these 

tables is chosen by the user with the restriction that all six coefficients 
must have the some dimension. 

For the sample cose pre.scntod, the coefficients are dimensioned 21 
X 11 and ore input in the following manner: 

"‘tI^ 

UtP ‘^Tl^ 

^ ^ ^ol5 * ^Tl^ 


*^T11^ 

9 

^^‘^al5’ ‘’‘TiP'* ‘’‘TlP 


7 
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III. PROGRAM OUTPUT 


A. Output Specifications 

The program allows the user to specify up to 90 vax’iables, con- 
tained in the first 999 locations of the common block, for printing or to 
be saved on tape. Printed output data arc assigned to Unit 6. Each 
block of printed data are preceded by the simulation time and case title. 
The block of printed data (6 variables per line of pxunt) uses an E15.8 
format with a four character identification name, 

All tape output data are written in binary form. Valuables saved 
for restart are assigned to Unit 13. Data to be plotted or used as input 
to other programs are assigned to Unit 12. The first variable of the 
data block saved is the simulation time in seconds. The remaining vari- 
ables, and their order of output, ax*o specified by the input data (See 
Section III, Code 2). 
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Appendix A 


PROGRAM VARIABtE DEFINITIONS 


The following program variables, located in the COMMON 

BLOCK , 

are alphabetized according to their Fortran Mnemonics. 


Program Symbol 

Variable Definition 


Units 

ABI(3,3) 

Transformation matrix from I 
frame 

to B- 

Unitless 

ACCOM(582) 

Dummy common 


Unitless 

ADG(3,3) 

Transibrmution matrix from G 
frame 

lo D- 

Unitless 

ADI(3,3) 

Transformation matrix from I 
frame 

to D~ 

Unitless 

ADR(3,3) 

Transformation matrix from R 
frame 

to D- 

Unitless 

AGR(3,3) 

Transformation matrix from R 
frame 

to G- 

Unitless 

AIB(3,3) 

Transformation matrix from B 
frame 

to I- 

Unitless 

AID(3,3) 

Transformation matrix from D 
frame 

to 1“ 

Unitless 

AIG(3,3) 

Transformation matrix from G 
frame 

to I" 

Unitless 

ALPHT 

Total angle-of-attack 


rad 

ARB(3,3) 

Transformation matx’ix iVom B 
frame 

to R- 

Unitless 

ARD(3,3) 

Transformation matrix irom D 
frame 

to R- 

Unities s 

ARI(3,3) 

Transformation matrix from I 
frame 

to R- 

Unities s 

ARP(3,3) 

Transformation matrix from P 
frame 

to R- 

Unitless 

ATCOM(45) 

Dummy common 


Unitless 

CA 

Axial force coefficient 


Unitless 

CCCOM(13) 

Dummy common 


Unitless 
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Program Symbol 

Voriuble Definition 

Units 

CH(3) 

Coefficients for geodetic latitude 
computations 

Unitless 

CEN 

Yawing moment coefficient 

Unitless 

CL 

Rolling moment coefficient 

Unitless 

CLMNT(3) 

Aerodynamic moment coefficients 
about vehicle center of mass 

Unitless 

CM 

Pitching moment coefficient 

Unitless 

CN 

Normal force coefficient 

Unitless 

CNV 

57.29577951 

deg/rad 

CR 

Oblate earth x'adial coefficient 

Unitless 

CW 

Oblate earth spin axis coefficient 

Unitless 

CPLMN(3) 

Aex'odynamic moment coefficient 
vector about moment reference 
point 

Unitless 

CXYZ(3) 

Acx’odynamic force coefficient 
vector 

Unitless 

CY 

Side force coefficient 

Unitless 

DIW(3) 

Dummy variable used in gravity 
gradient tox’quo calculation 

1 

slugs- ft' 

DPSl 

Difference between geodetic and 
geocentxac latitudes 

rad 

DREF 

Aex’odynamic reference diameter 

ft 

DT 

Integration time step 

sec 

DTP 

Output frequency for tape save 

sec 

DTSAM 

Specified time to save vaxiablcs 
for restart 

sec 

EPCOM(66) 

Dummy common 

Unitless 

FAB (3) 

B-frame acrody?vamic tbx’ces 

lb 

FAI(3) 

I-frame aerodynamic fox’ces 

lb 

FAR(3) 

ll-frame aerodynamic forces 

lb 

FARC (3) 

R-frame aex’odynamic force coeffi- 
cients 

Unitless 

FCA(21,11) 

Table of axial force coefficients 

Unitless 

FCEN(21,11) 

Table of yawing moment coeffi- 
cients about moment x’cference 
point 

Unitless 
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Proifrom Symbol 

Vaylablo Doitnitlon 

Units 

FCL(21,U) 

Table of i*olling moment coefficients 
about moment reference point 

Unitless 

FCM(21,11) 

Table of pitching moment coofflci*’ 
ents about moment reference point 

Unitloss 

FCN(21,11) 

Table of normal force cooffleionts 

Unitlcss 

FCY(21,11) 

Table of side force coefficients 

Unitless 

FOMK3) 

Aerodynamic term in translational 
acceleration calculation 

ft/sec^ 

FRHO(ll,37) 

Table of atmospheric densities 

kg/m^ 

G(3) 

Gravitational acceleration 

ft/soc^ 

GAMI 

Inertial flight path angle (positive 
up fron) local geocentric horizontal) 

rod 

GAMR 

Relative flight path angle (positive 
up from local geodetic horizontal) 

rad 

GCR 

Temporary variable in acceleration 
of gravity calculation 

ft/sec^ 

GGB(3) 

Gravity gradient toi’ques in 
frame 

ft-lb 

GMAG 

Temporary variable in acceleration 
of gravity calculation 

ft/sec^ 

GMASI 

Reciprocal of vcliicle mass 

1 /slugs 

GMASS 

Veliicle mass 

slugs 

H 

Geodetic altitude 

ft 

HB(3) 

Angular momentum in B -frame 

ft“lb~sec 

IXYZ(3,3) 

Moment of inertia tensor 

2 

slugs “ft 

IXYZI(3,3) 

Inverse of moment of inertia tensor 

l/slugS“ft 

LAMDA 

Inertial longitude measured in I- 
frame 

I’ad 

LAMDE 

Earth fixed loiigitude 

rad 

LMRP(3) 

Temporary vaxiable in aerodynamic 
moment calculation 

Unitless 

OMGE 

Earth's I’otation rate 

rad /sec 

PCDUM(12) 

Dummy common 

Unitless 

PHIA 

Aci’odynamic roll angle 

rad 
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Proernm Symbol 

Variable Definition 

Units 

PHIBK 

Bank angle 

rad 

PI 

tr 

rud 

^qiU3) 

Angular rates in B-framo 

rad /sec 

FQRD(3) 

Angular accelerations in B “frame 

rad/sec^ 

PSl 

Geocentric latitude (positive north, 
negative south of equator) 

rad 

PS ID 

Geodetic lotitudc (positive north, 
negative south of equator) 

rad 

Q(4) 

Quaternion parameters 

Unitless 

QBAR 

Dynamic pressure 

Ib/ft^ 

QD(4) 

Time derivities of quaternions 

1/sec 

QSQR 

/ Q(l)^ + Q(2)^ + Q(3)^ Q(4)^ 

Unitless 

R(3) 

Vehicle I -frame position 

ft 

RX(3) 

Unit vector along vehicle I-frame 
position 

Unitless 

RB1(3) 

Unit vector along vehicle B -frame 
position 

Unitless 

RDOT(3) 

Vehicle I-framc translational velocity 

ft /sec 

RHO 

Atmospheric density 

slugs /ft^ 

RMAG 

Magnitude of vehicle position vector 

ft 

RPSI 

Radius of the earth 

ft 

RR,RR2 

Temporary variables used in 
gravitational acceleration 
calculations 

Unitless 

SD.SJ^SH 

Earth's gravitational constants 

Unitless 

SIGI 

Inertial heading (positive clockwise 
from north) 

rad 

SIGR 

Relative heading (positive clockwise 
from north) 

rad 

SPSI,S?SI2 

Temporary variables used in gravita- 
tional acceleration calculation 

Unitless 

SREF 

Aerodynamic reference area 

ft^ 

SWD,SWH 

Temporary variables used in gravi- 

Unitless 


tational acceleration calculation 
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Procrram Symbol 

Vai'inble Deflnilion 

Units 

TAB (3) 

Total oxtoi'nal torques in B- frame 

ft-lb 

TABP(U) 

Table of anglo-of-attack values 

rad 

TEND 

Simulation stop time 

see 

THOR 

Simulation time 

hr 

TIME 

Simulation time 

sec 

TLAT(U) 

Table of geocentric latitude values 

rad 

TLNG(37) 

Table of longitude values 

rod 

TOA(3) 

Aerodynamic torques in B-frame 

ft»lb 

TPHIA(21) 

Table of aerodynamic roll angle 
values 

rad 

TRUN 

Total run time 

sec 

TSUM(3) 

Time rate of ebange of vehicle 
angular momentum 

ft lb 

TWOPI 

2 !t 

rad 

V(3) 

Inertial velocity in I“frnme 

ft/sec 

VACOM<67) 

Dummy common 

Unitless 

VATM(3) 

Atmospheric velocity in I-frame 

ft /sec 

VCCOM(18) 

Dummy common 

Unitless 

VDCOM(X07) 

Dunuay common 

Unitless 

VD'0T(3) 

Inertial translation acceleration in 
I “frame 

ft/sec^ 

VID(3) 

Inertial velocity in D-* frame 

ft /sec 

VIGO) 

Inertial velocity in G-frame 

ft /sec 

VIMAG 

Magnitude of inertial velocity 

ft /sec 

VPCOM(39) 

Dummy common 

Unitless 

VRELBO) 

Relative velocity in B -frame 

ft /see 

VRELDO) 

Relative velocity in D-frarne 

ft /sec 

VRELGO) 

Relative velocity in G-frame 

ft /sec 

VRELI(3) 

Relative velocity in 1- frame 

ft /sec 

VRMAG 

Magnitude of relative velocity 

ft /sec 

WT 

Vehicle weight 

lb 

WXH(3) 

Angular momentum direction change 

ft-lb 

XCG,YCG,ZCG 

Cg location in N-frame along x,y, 
and z-axes, respectively 

ft 


Program Symbol 

XDUMC27) 

XDDUMC27) 

XMRP, YMRP, ZMRP 


Variable Definition 


Units 


Dummy common Unitless 

Dummy common Unitless 

Moment reference point in N-* frame ft 

along X, y, and z-axes, 

respectively 
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Appendix B 


REFERENCE FRAMES AND TRANSFORMATIONS 


A . Reference Frames 

All reference frames are right handed systems (Figs. B-1, B-2, 
and B-3). 


1. I — Inertial Frame . The I-frame has its origin at the center 
of the earth with the Xj axis through the Greenwich mex’idian at time 

zero. The Zj axis points through the North Pole and the Yj axis com- 
pletes the right handed system. It is in this frame the accelerations are 
integrated . 

2. G — Geocentric Frame. The G-frame has its origin at the 
vehicle's center of mass. The Xq axis points north, the axis points 

east, and the axis is directed downward along the radius vector to 

the earth's center. 

3. D — Geodetic Frame. The D-frame has its origin at the vehicle's 
center of mass (Fig. B-2). The Xj^ axis points north, the Yjj axis points 

east, and the Z^ axis is directed downward along the local geodetic. 

4. R — Relative Velocity Frame . The R-frame has its origin at 
the vehicle's center of mass (Fig. B-2”). The Xj^ axis is directed along 

the relative velocity vector. The axis is directed downward in a 

plane containing the velocity vector and the local geodetic. The Yj^ 

axis completes the right handed system. 

5. B — Body Frame. The body fixed B -frame has its origin at the 
vehicle's center of mass (Fig. B-3). The direction of the axes are chosen 
so as to be consistent with the definition of aerodynamic parameters. It 
is in this frame that the external forces and moments are computed. 

6. N — Input Data Reference Frame. This body fixed reference 
frame is parallel to the B-frame coordinate system with its origin chosen 
by the user. It is in this frame that the aerodynamic data and mass 
properties data are read into the program. Usually, it is convenient to 
choose the aerodynamic moment reference point as the origin of this frame. 
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EQUATORIAL PLANE 



A®'-[-(90+^^)J2 I M 3 

Figure B-1. Inertial (I), geocentric (G) and geodetic (D) 

reference frames. 



Jodetic (D) and relative velocity (H) 
reference frames. 



B . Transformations 


The transformation matrix transforms vectors from the L-frame 
to the M-frame. The matrices are formed by successive rotations through 
the indicated Euler angles. The rotation [6]. is used to indicate the 

direction cosine matrix for a positive rototion about the i^" axis through 
the angle 6. The sequence of rotations is read from right to left. Since 
all transformations shown are orthogonal, the inverse transformation 
matrices are merely the transpose of those given. 

C* I 

1. Transformation from I to G-frame A 

A®* = l-(90 + [Alg 

where X = inertial longitude 
\p = geocentric latitude 

2. Transformation from I to D-frame A^^ 

A^* = [-(90 + i|ij3)]2 [Xlg 

where X = inertial longitude 
= geodetic latitude 

3. Transformation from D to R-frame A^^ 


A"*' = fYRl2[aRl3 

where = relative heading 

= relative flight path angle 


4. Transformation from R to B-frame 


where (p 


3 


a„ 


a 


= bank angle 

= total angle-of-attack 
= aerodynamic roll angle 


5. Transformation from I to B -frame 

Initially, the transformation from the I to B-frame is computed 

by A»“ a“‘ 

From the inverse of this transformation, the initial quaternion 
parameters are computed. The four quaternion parameters are defined 

as follows: 

Qj = a sin (<|i/2) 


Qg = P sin ((|>/2) 


Qj s Y sin (<l'/2) 


= cos ((t>/2) 


where o, 3, and y eigenaxis direction cosines, and <p is the eigenaxis 
rotation angle. The quaternions are initialized as follows: 




®11 ®12 ®13 


“21 ®22 ®23 


®31 ®32 ®33 


Q 4 = j / 1 + a^^ + 022 ^33 


KQ4 i 0, 


Qi = 


(^QO "" 


1 " '“32 " “ 23 ^ 


Qg = (ai 3 - Ogi) 

^3 “ ^“21 “ “12^ 

If Q4 = 0 , 

Ql = >/ - j (“22 ■*■ ®33^ 

If 0, 

^2 ~ ^^12 ■*■ *^2p 

*^3 " ^“i3 "31^ 

If Ql = 0. 

Q 2 = 4 (^22^ 

If Q2 ^ 0 , 

^3 “ 4^ ^“23 ^32^ 

If Qg = 0, 

Q3 = 1 

From the integrated values of the quaternion parameter time deriva- 
IB 

ties, the A transformation matrix is updated as shown in Appendix C. 
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Appendix C 


EQUATIONS OF MOTION 


This appendix presents briefly the equations of motion used to 
describe the dynamics of a vehicle in 6-degrees-of-freedom. The equa- 
tions were derived for a vehicle of constant moss distribution and the 
external forces and torques considered were those resulting only from 
gravity end aerodynamics. The mathematical models for the Earth and 
gpravity can be found in Reference 1. The atmospheric density is an 
input to the program as a function of position in orbit. 

The aerodynamic force and moment coefficients and the atmospheric 
density are determined by table lookup at each integration step. The 
aerodynamic coefficients (C^, Cy, C^) are a function of 

total angle-of-attack, a,j,, and aerodynamic roll angle, The density 

(p) is a function of geocentric latitude, and inertial longitude, 

These angles are defined as 


a,j, = tan 


-1 


A 


'^RELB^^^ 


0 < a,p 


< 180° 


. - tnn"^ '^RELB^^^ 


RELB' 


0 < 4) < 360° 

— — 


ij; = tan 


-1 


11(3) 


/ R(l)^ + R(2)^ 


- 90° < ijj < 90° 


X = tan 


-1 R(2) 
R(T) 


0 < A < 360° 


where 


VjjElb “ relative velocity vector in B- frame 


,R = inertial I- frame position vector. 
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A . Translational Motion 

The translationa) equation of motion for a body in the inertial 
I~foame is 


R 




where 



V = inertial acceleration vector 
gravitational acceleration vector 
mass of the body 
= I-frome aerodynamic forces 


The I-frame velocity equation is 


n = V ^ R^ + /Rdt 


where 


^o “ X) “ initial I-framo velocity vector. 

Integrating the velocity vector and adding the initial position, the 
I-frame position vector is 

R = R^ + / R dt 

where 


initial I-frame position veetor. 


The aerodynamic forces are computed in the body B -frame and 
tranaformed into the I-frame for the translational acceleration computation. 
The B'frame aerodynamic forces are 


^AB “ *1bAR ®Ref ^YZ 


where 


'XYZ 


**BAR ” 


-C. 


-C 


N 


= B-frame aerodynamic force 
coefficient vector 


S 


Ref 


^ P 2 

aerodynamic reference area. 


RELB I s dynamic pressure 


Then the 1-frame aerodynamic forces are 



A*® F 


AB 


B. Rotational Motion 

The equation describing the rotational motion about the center of 
mass of a rigid body in the rotating B-frame is given by 


[I] 0) = T^g - w X ([I] .w) 

where 



*xx 

■*XY 

■*xz 


[I] = 

“^XY 

lyY 

“*YZ 

= inertia tensor 


“*XZ 

"*YZ 

^ZZ 
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w = PQRD ss B~frame ong;ular aceleration vector 
w s PQIt ss B-frame angular velocity vector 
*^AB s B-frame external torques. 

Then 


w = + / tU"^ tT^g - w X ([I].w) dt 

where 

Wq - initial B-frame angular velocity vector. 


The aerodynamic torques about the vehicles center of mass are 


^OA “ '^BAR ®Ref ^Ref 


^Li 


MN 


^^MRP ■* 

SyzJ 


where 


D 


Ref 


aerodynamic reference diameter 


"PLMN = 


m 


n 


moment coefficients about the moment 
reference point 


Xjy^Rp - aerodynamic moment reference point 
^CG “ 


center of mass location. 


The gravitational torques are given by 


tRJb - tHffl X [II RBi] 

1B|® 


where 


GM s product of universal gravitational constant and the mass 
of the earth 

RFi = unit vector along body B-frame position. 

The B-frame external torques are 

*AB = *OA * SBff , 

The quaternion parameters are updated in a simular manner: 

Q = + / q dt 


where 


= initial quaternion parameters defined in Appendix B. 
The time derivity Q is defined as 


1 

1 

II 


m 

■**3 ^2 




Q4-Q1 

^3 


-«2 

Q1-Q4 




"'®2 "^3 
■ 
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From these integrated values of the A*® transformation matrix 
is calculated as follows: 

(2 q/ + 2 2(Qj^ Q2 - Qg Q4) 2(Qj^ Qg + Qg Q^) 

2 (Qi Qg + Qg Q4) (2 + 2 q/ - 1 ) 2 (Qg Qg - Q4) 

2(Qj Qg - Qg Q4) 2(Q2 Qg + Qj Q^) (2 Qg^ + 2 Q^^ -1) 
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